clear all
close all
clc

% load data: red, green, nlue, cyan, magenta, yellow, rgmix

load('filter_data.mat');
load('D65.mat');
load('color_matching_func.mat');

% set D65 params
d65 = d65source(1:10:end,2);

% set color matching function params
cmf = color_matching_func(1:10:391,:);

xbar = cmf(:,2);
ybar = cmf(:,3);
zbar = cmf(:,4);

%% Plot D65 Spectral Power Distribution

figure

hold on
plot(d65source(:,1),d65source(:,2),'k','LineWidth',2);
title('Spectral Power Distribution of CIE D65')
xlabel('Wavelength (nm)')
ylabel('Power')
hold off

print('-depsc', 'd65_spec.eps')

%% Plot Transmission Spectrum

figure

hold on
plot(red(:,1),red(:,2),'r','LineWidth',2);
plot(green(:,1),green(:,2),'g','LineWidth',2);
plot(blue(:,1),blue(:,2),'b','LineWidth',2);
plot(cyan(:,1),cyan(:,2),'c','LineWidth',2);
plot(magenta(:,1),magenta(:,2),'m','LineWidth',2);
plot(yellow(:,1),yellow(:,2),'y','LineWidth',2);
plot(rgmix(:,1),rgmix(:,2),'k','LineWidth',2);
title('Spectral Transmittance of Filters')
xlabel('Wavelength (nm)')
ylabel('Percent Transmission')
legend('red','green','blue','cyan','magenta','yellow','red-green mix','Location','EastOutside')
hold off

% saved as 'filter_transmittance.eps'

%% Calculate xy Chromaticity Coordinates

k = 100/sum(d65.*ybar);

% red filter xy
X_red = k*sum(xbar.*(red(:,2)/100).*d65);
Y_red = k*sum(ybar.*(red(:,2)/100).*d65);
Z_red = k*sum(zbar.*(red(:,2)/100).*d65);

x_red = X_red/(X_red + Y_red + Z_red);
y_red = Y_red/(X_red + Y_red + Z_red);

% green filter xy
X_green = k*sum(xbar.*(green(:,2)/100).*d65);
Y_green = k*sum(ybar.*(green(:,2)/100).*d65);
Z_green = k*sum(zbar.*(green(:,2)/100).*d65);

x_green = X_green/(X_green + Y_green + Z_green);
y_green = Y_green/(X_green + Y_green + Z_green);

% blue filter xy
X_blue = k*sum(xbar.*(blue(:,2)/100).*d65);
Y_blue = k*sum(ybar.*(blue(:,2)/100).*d65);
Z_blue = k*sum(zbar.*(blue(:,2)/100).*d65);

x_blue = X_blue/(X_blue + Y_blue + Z_blue);
y_blue = Y_blue/(X_blue + Y_blue + Z_blue);

% cyan filter xy
X_cyan = k*sum(xbar.*(cyan(:,2)/100).*d65);
Y_cyan = k*sum(ybar.*(cyan(:,2)/100).*d65);
Z_cyan = k*sum(zbar.*(cyan(:,2)/100).*d65);

x_cyan = X_cyan/(X_cyan + Y_cyan + Z_cyan);
y_cyan = Y_cyan/(X_cyan + Y_cyan + Z_cyan);

% magenta filter xy
X_magenta = k*sum(xbar.*(magenta(:,2)/100).*d65);
Y_magenta = k*sum(ybar.*(magenta(:,2)/100).*d65);
Z_magenta = k*sum(zbar.*(magenta(:,2)/100).*d65);

x_magenta = X_magenta/(X_magenta + Y_magenta + Z_magenta);
y_magenta = Y_magenta/(X_magenta + Y_magenta + Z_magenta);

% yellow filter xy
X_yellow = k*sum(xbar.*(yellow(:,2)/100).*d65);
Y_yellow = k*sum(ybar.*(yellow(:,2)/100).*d65);
Z_yellow = k*sum(zbar.*(yellow(:,2)/100).*d65);

x_yellow = X_yellow/(X_yellow + Y_yellow + Z_yellow);
y_yellow = Y_yellow/(X_yellow + Y_yellow + Z_yellow);

% d65 xy

X_d65 = sum(xbar.*d65);
Y_d65 = sum(ybar.*d65);
Z_d65 = sum(zbar.*d65);

x_d65 = X_d65/(X_d65 + Y_d65 + Z_d65);
y_d65 = Y_d65/(X_d65 + Y_d65 + Z_d65);

% Spectrum Locus

xlocus = xbar./(xbar + ybar + zbar);
ylocus = ybar./(xbar + ybar + zbar);

xlocus_smooth = color_matching_func(:,2)./(color_matching_func(:,2) + color_matching_func(:,3) + color_matching_func(:,4));
ylocus_smooth = color_matching_func(:,3)./(color_matching_func(:,2) + color_matching_func(:,3) + color_matching_func(:,4));

%% Plot Spectral Locus and xy coordinates

figure

labels = cellstr( num2str(cmf(:,1)) );
labelsleft = cellstr( num2str(cmf(14:17,1)) );
labelsright = cellstr( num2str(cmf(18:26,1)) );

hold on
% plot locus
plot(xlocus, ylocus, '.k', 'MarkerSize', 10);
plot(xlocus_smooth, ylocus_smooth, '-k', 'LineWidth', 2);
plot([xlocus(1),xlocus(end)], [ylocus(1), ylocus(end)], 'k', 'LineWidth', 2);
%plot points
plot(x_red, y_red, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
plot(x_green, y_green, 'o', 'MarkerFaceColor', 'g', 'MarkerSize', 8);
plot(x_blue, y_blue, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 8);
plot(x_cyan, y_cyan, 'o', 'MarkerFaceColor', 'c', 'MarkerSize', 8);
plot(x_magenta, y_magenta, 'o', 'MarkerFaceColor', 'm', 'MarkerSize', 8);
plot(x_yellow, y_yellow, 'o', 'MarkerFaceColor', 'y', 'MarkerSize', 8);
plot(x_d65, y_d65, 'o', 'MarkerFaceColor', 'k', 'MarkerSize', 8);

% text(xlocus(7), ylocus(7), labels{7}, 'VerticalAlignment','top','HorizontalAlignment','right')
% text(xlocus(12), ylocus(12), labels{12}, 'VerticalAlignment','top','HorizontalAlignment','right')
% text(xlocus(13), ylocus(13), labels{13}, 'VerticalAlignment','top','HorizontalAlignment','right')
% text(xlocus(14:17), ylocus(14:17), labelsleft, 'VerticalAlignment','bottom','HorizontalAlignment','right')
% text(xlocus(18:26), ylocus(18:26), labelsright, 'VerticalAlignment','bottom','HorizontalAlignment','left')
% text(xlocus(28), ylocus(28), labels{28}, 'VerticalAlignment','bottom','HorizontalAlignment','left')
% text(xlocus(end), ylocus(end), labels{end}, 'VerticalAlignment','top','HorizontalAlignment','left')

title('Chromaticity Plot for Filters');
xlabel('x');
ylabel('y');
axis([-0.1 0.9 -0.1 0.9]);
set(gca,'YTick',[-0.1, 0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
set(gca,'XTick',[-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
box on
grid on
axis square
hold off

print('-depsc', 'filter_xy.eps')

%% Determining the dominant wavelength HERE%%%%%%%%%%%%%%%%%%%%%%%%%%
% just need the x chromaticity coordinate for the dominate wavelength for
% each color: x_dom_wave_red, x_dom_wave_green, x_dom_wave_blue, etc...

slope_red = (y_red - y_d65)/(x_red - y_d65);
yint_red = y_red - slope_red*x_red;

slope_green = (y_green - y_d65)/(x_green - y_d65);
yint_green = y_green - slope_green*x_green;

slope_blue = (y_d65 - y_blue)/(y_d65 - x_blue);
yint_blue = y_blue - slope_blue*x_blue;

slope_cyan = (y_d65 - y_cyan)/(y_d65 - x_cyan);
yint_cyan = y_cyan - slope_cyan*x_cyan;

slope_magenta = (y_magenta - y_d65)/(x_magenta - y_d65);
yint_magenta = y_magenta - slope_magenta*x_magenta;

slope_yellow = (y_yellow - y_d65)/(x_yellow - y_d65);
yint_yellow = y_yellow - slope_yellow*x_yellow;

figure

hold on
% plot locus
plot(xlocus, ylocus, '.k', 'MarkerSize', 10);
plot(xlocus_smooth, ylocus_smooth, '-k', 'LineWidth', 2);
plot([xlocus(1),xlocus(end)], [ylocus(1), ylocus(end)], 'k', 'LineWidth', 2);
% plot lines for dominant wavelength
plot([x_d65, 0.8], [y_d65,slope_red*0.8 + yint_red])
plot([x_d65, 0.5], [y_d65,slope_green*0.5 + yint_green])
plot([x_d65, 0], [y_d65,slope_blue*0 + yint_blue])
plot([x_d65, 0], [y_d65,slope_cyan*0 + yint_cyan])
plot([0.25,x_d65, 0.4], [slope_magenta*0.25 + yint_magenta,y_d65,slope_magenta*0.4 + yint_magenta])
plot([x_d65, 0.5], [y_d65,slope_yellow*0.5 + yint_yellow])
% plot points
plot(x_red, y_red, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
plot(x_green, y_green, 'o', 'MarkerFaceColor', 'g', 'MarkerSize', 8);
plot(x_blue, y_blue, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 8);
plot(x_cyan, y_cyan, 'o', 'MarkerFaceColor', 'c', 'MarkerSize', 8);
plot(x_magenta, y_magenta, 'o', 'MarkerFaceColor', 'm', 'MarkerSize', 8);
plot(x_yellow, y_yellow, 'o', 'MarkerFaceColor', 'y', 'MarkerSize', 8);
plot(x_d65, y_d65, 'o', 'MarkerFaceColor', 'k', 'MarkerSize', 8);

% text(xlocus(7), ylocus(7), labels{7}, 'VerticalAlignment','top','HorizontalAlignment','right')
% text(xlocus(12), ylocus(12), labels{12}, 'VerticalAlignment','top','HorizontalAlignment','right')
% text(xlocus(13), ylocus(13), labels{13}, 'VerticalAlignment','top','HorizontalAlignment','right')
% text(xlocus(14:17), ylocus(14:17), labelsleft, 'VerticalAlignment','bottom','HorizontalAlignment','right')
% text(xlocus(18:26), ylocus(18:26), labelsright, 'VerticalAlignment','bottom','HorizontalAlignment','left')
% text(xlocus(28), ylocus(28), labels{28}, 'VerticalAlignment','bottom','HorizontalAlignment','left')
% text(xlocus(end), ylocus(end), labels{end}, 'VerticalAlignment','top','HorizontalAlignment','left')

title('Chromaticity Plot for Filters With Dominant Wavelength Line');
xlabel('x');
ylabel('y');
axis([-0.1 0.9 -0.1 0.9]);
set(gca,'YTick',[-0.1, 0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
set(gca,'XTick',[-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
box on
grid on
axis square


hold off

print('-depsc', 'filter_dom_wave_xy.eps')

%% Calculating the Purity

x_dom_wave_red = 0.6826;
x_dom_wave_green = 0.3445;
x_dom_wave_blue = 0.08708;
x_dom_wave_cyan = 0.5932;
x_dom_wave_magenta = 0.2873;
x_dom_wave_yellow = 0.4441;

% red
purity_red = (x_red - x_d65)/(x_dom_wave_red - x_d65);

% green
purity_green = (x_green - x_d65)/(x_dom_wave_green - x_d65);

% blue
purity_blue = (x_blue - x_d65)/(x_dom_wave_blue - x_d65);

% cyan
purity_cyan = (x_cyan - x_d65)/(x_dom_wave_cyan - x_d65);

% magenta
purity_magenta = (x_magenta - x_d65)/(x_dom_wave_magenta - x_d65);

% yellow
purity_yellow = (x_yellow - x_d65)/(x_dom_wave_yellow - x_d65);

%% Additively Mixing Tristimulus Values for xy coordinates

% new tristimulus values

X_new_rgmix = (X_red + X_green)/2;
Y_new_rgmix = (Y_red + Y_green)/2;
Z_new_rgmix = (Z_red + Z_green)/2;

% xy coordinates

x_new_rgmix = X_new_rgmix/(X_new_rgmix + Y_new_rgmix + Z_new_rgmix);
y_new_rgmix = Y_new_rgmix/(X_new_rgmix + Y_new_rgmix + Z_new_rgmix);

% dominant wavelength

slope_new_rgmix = (y_new_rgmix - y_d65)/(x_new_rgmix - y_d65);
yint_new_rgmix = y_new_rgmix - slope_new_rgmix*x_new_rgmix;

%% Subtractively "mixing" Spectra for xy coordinates

rgmix_spectra = red(:,2).*green(:,2);

% xy coordinates
% red-green spectrally computed overlapped filters xy
X_rgmix_spectra = k*sum(xbar.*(rgmix_spectra/100).*d65);
Y_rgmix_spectra = k*sum(ybar.*(rgmix_spectra/100).*d65);
Z_rgmix_spectra = k*sum(zbar.*(rgmix_spectra/100).*d65);

x_rgmix_spectra = X_rgmix_spectra/(X_rgmix_spectra + Y_rgmix_spectra + Z_rgmix_spectra);
y_rgmix_spectra = Y_rgmix_spectra/(X_rgmix_spectra + Y_rgmix_spectra + Z_rgmix_spectra);


% dominant wavelength

slope_rgmix_spectra = (y_rgmix_spectra - y_d65)/(x_rgmix_spectra - y_d65);
yint_rgmix_spectra = y_rgmix_spectra - slope_rgmix_spectra*x_rgmix_spectra;

%% Measurement of superimposed filters

% xy coordinates
% red-green overlapped filters xy
X_rgmix = k*sum(xbar.*(rgmix(:,2)/100).*d65);
Y_rgmix = k*sum(ybar.*(rgmix(:,2)/100).*d65);
Z_rgmix = k*sum(zbar.*(rgmix(:,2)/100).*d65);

x_rgmix = X_rgmix/(X_rgmix + Y_rgmix + Z_rgmix);
y_rgmix = Y_rgmix/(X_rgmix + Y_rgmix + Z_rgmix);

% dominant wavelength

slope_rgmix = (y_rgmix - y_d65)/(x_rgmix - y_d65);
yint_rgmix = y_rgmix - slope_rgmix*x_rgmix;

% chromaticities
rg_xy = [x_rgmix,y_rgmix;x_rgmix_spectra,y_rgmix_spectra;x_new_rgmix,y_new_rgmix];
% make curves!

%% Plot for mixtures

figure

hold on
% dummy points for legend
plot(-3,-3, 'o', 'MarkerEdgeColor','k','MarkerFaceColor', [0.7 0.5 0], 'MarkerSize', 8);
plot(-3,-3, 'd','MarkerEdgeColor','k', 'MarkerFaceColor', [0.7 0.5 0], 'MarkerSize', 8);
plot(-3,-3, 's','MarkerEdgeColor','k', 'MarkerFaceColor', [0.7 0.5 0], 'MarkerSize', 8);
% locus
plot(xlocus, ylocus, '.k', 'MarkerSize', 10);
plot(xlocus_smooth, ylocus_smooth, '-k', 'LineWidth', 2);
plot([xlocus(1),xlocus(end)], [ylocus(1), ylocus(end)], 'k', 'LineWidth', 2);
% dominant wavelength
plot([x_d65, 0.8], [y_d65,slope_new_rgmix*0.8 + yint_new_rgmix])
plot([x_d65, 0.8], [y_d65,slope_rgmix_spectra*0.8 + yint_rgmix_spectra])
plot([x_d65, 0.8], [y_d65,slope_rgmix*0.8 + yint_rgmix])
% xy coordinates of mixtures
plot(x_rgmix,y_rgmix, 'o', 'MarkerEdgeColor','k','MarkerFaceColor', [0.7 0.5 0], 'MarkerSize', 8);
plot(x_rgmix_spectra,y_rgmix_spectra, 'd','MarkerEdgeColor','k', 'MarkerFaceColor', [0.7 0.5 0], 'MarkerSize', 8);
plot(x_new_rgmix,y_new_rgmix, 's','MarkerEdgeColor','k', 'MarkerFaceColor', [0.7 0.5 0], 'MarkerSize', 8);
% xy coordinates of original filters
plot(x_red, y_red, 'o', 'MarkerFaceColor', 'r', 'MarkerSize', 8);
plot(x_green, y_green, 'o', 'MarkerFaceColor', 'g', 'MarkerSize', 8);
% xy coordinates of source
plot(x_d65, y_d65, 'o', 'MarkerFaceColor', 'k', 'MarkerSize', 8);
title('Chromaticity Plot for Mixtures With Dominant Wavelength Line');
xlabel('x');
ylabel('y');
axis([-0.1 0.9 -0.1 0.9]);
set(gca,'YTick',[-0.1, 0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
set(gca,'XTick',[-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9])
legend('Measured','Subtractive','Additive','Location','NorthEast')
box on
grid on
axis square

hold off

print('-depsc', 'mixture_dom_wave_xy.eps')

% purity for mixtures

x_dom_wave_new_rgmix = 0.581;
x_dom_wave_rgmix_spectra = 0.613;
x_dom_wave_rgmix = 0.6402;

purity_new_rgmix = (x_new_rgmix - x_d65)/(x_dom_wave_new_rgmix - x_d65);
purity_rgmix_spectra = (x_rgmix_spectra - x_d65)/(x_dom_wave_rgmix_spectra - x_d65);
purity_rgmix = (x_rgmix - x_d65)/(x_dom_wave_rgmix - x_d65);

